Last month, I was reminded at work that I wanted to write a fluid simulator. In 2018, just before Christmas, I was allowed to spend a week at work to try to make fluid simulation work, since we were in-between projects. A guy wrote an algorithm in 2D for a game my company did, using a chapter coming from GPU Gems books. It might have been this one, referencing Robert Bridson's book: Real-time simulation and rendering of 3D fluids. This one feels similar too: Fast fluid dynamics simulation on the GPU
The guy told me it was based on Robert Bridson's book, so I tend to believe it was the first of the two. I decided to buy the book in 2018 : Fluid Simulation for Compute Graphics.
There are other implementations that one could find for free online, but I wanted to write my own. The idea is to solve a simplified version of Navier-Stokes equations:
$$ { \begin{aligned} \frac{\partial {\vec{u}}}{\partial {t}} + \vec{u} \cdot \nabla{\vec{u}} + \frac{1}{\rho} \nabla{p} = \nu\nabla \cdot \nabla{\vec{u}} \\ \nabla \cdot \vec{u} = 0 \end{aligned} } $$
The book contains all the information to solve these equations within the first five chapters, ending at page 97. I did read the chapter 4, but did not program it yet since it is about collisions with the environment, which I do not take into account yet. This makes for an incredibly dense book, where it might be necessary to read a bit more outside of it when needed.
The book uses the first two chapters to simplify the problem and split it into four steps. Only two of them are complex, and only one of those two needs a linear solver. They are listed in Chapter 2 of the book at page 20. Each of them is repeated for each frame.
Determining a good time step is just there to ensure that the frame duration is short enough for the fluid velocity, otherwise it would move too fast compared to the size of each cells. It is mostly a rule of thumb and easy to compute. Applying external forces is also very easy. Just apply gravity on each velocity using \( \vec{u_{n+1}} = \vec{u_{n}} + \Delta t \space \vec{g} \).
The two real steps are the advection and the pressure solve (the pressure solve being way harder). I started on that in 2018. I decided to use the Octave programming language (can we call it that?) since I had prior experience with MATLAB, and it is a free clone of it. It has the linear solver libraries for when I need it, and visualisation tools too.
Before even talking about the advection and pressure, I need to talk a bit about interpolation. I thought it would be an easy step, but nope, it was way harder than expected. The reason is the way the data is stored, using special grids names MAC grids (Marker-and-cell, I don't understand why it is called like that). The idea is to store most values (most importantly the pressure) in the center of the cells, but the velocites on the edges of the cells.
The matrix holding the pressure does not have the same dimensions as the one holding the X velocities, which in turn does not have the same dimensions as the one holding the Y velocities. Thus, it is very hard to keep track of the indices properly, let alone the interpolations. The same exact algorithm cannot even be used to interpolate the pressure and the velocities themselves (what a pain). This is why I spent some time testing the interpolation.
For now, I only use linear interpolation (bilinear in 2D) meaning that if a value is 20% between A and B (closest to A), it will have 80% of the value in A and 20% of the value of B. For better results, Bridson mentions we should use cubic interpolation (bicubic in 2D or tricubic in 3D). I probably could use other stuff like Hermite or cubic splines, but for now, I am good enough with linear. Here is a picture I found on Wikipedia ( Bicubic interpolation ) :
On grids having values on different places on the grid, tricubic interpolation becomes very slow, but we have to live with that to have a good precision. Velocity interpolation is especially annoying because the X and Y components of the velocities are not at the same locations on the grid. Just showing how much of an impact bicubic interpolation can have over bilinear:
I will stop now, but I will keep that in mind for the future. Every time I query a value in the grid at any position, I have to interpolate the value like this, so it is a critical step.
The first of the main steps is the advection, covered by chapter 3 of Bridson's book. This is conceptually simple: the idea is to move the physical values in the direction of the velocities, just like we would do for particles. Heat will move with the hot hair, just like pressure and smell for example. The trick here, is that we want pressure and velocities at specific locations in the grid. So instead of moving particles forward, we suppose the grid points are the forward points and we try to find where that particle was last frame by moving in the reverse direction of the particle. Once we find (thanks to velocity) where it was last frame, we compute (thanks to interpolation) what the values were at that location last frame and set them at the grid location we are working on. It's just logical that pressure (and heat, and velocity itself) move with the current (or the wind).
To illustrate the algorithm, I found this picture in one the GPU gems papers mentioned earlier.
Obviously, the velocity in our case is not stored in the center of the cells, but we get the point. The picture even shows two more things. First, from which points to interpolate (in the pressure case, it's like that, but not for velocities themselves). Second, the equation used to find where the particle was last frame : \( -\vec{u_n} \space \Delta t \).
The equation is actually a numerical integration, called Euler integration (in this case, it is a forward Euler integration even though we go backwards).
There are more accurate integration algorithms, the most well-known being Runge-Kutta. Different versions of the algorithm exist (order 2, order 4). In those algorithms, the idea is to integrate somewhere in the middle, compute the velocity there, and then use those middle velocities to integrate, since they are closer to the actual velocity of the particle during the time step. For more details : Wikipedia on Runge-Kutta . In conclusion, I will probably have to go implement this properly in the future, but for now I do the simple thing. See the next picture to understand how far we can get from the real values by using the simplest Euler equation.
I just implemented that logic and ran it in a loop and looked at the results for different values. The velocity advection by itself is interesting, like shown just below. We can see that the velocity is pointing outwards and advects itself out of the simulation.
Now for the hard part of the solve. This is where I failed in 2018. I had something that was kind of working, but not really at the same time... and I ended up failing to do it all in 4 days (pretty sure it was a 4 day week, just before Christmas).
The code I had at the time was saved somewhere on a Perforce repository, and it was probably deleted not long after. Obviously I could not save the code for personal use, so I just lost it, in a very incomplete state. Since them, I kept thinking about starting over, but never actually did it. I already planned around 2022 or 2023 to incorporate some form of fluid simulation in my engine. At work, another employee asked me lately in front of the team's intern if I was missing doing real physics, as I did before. I answered that I could still work on real physics, and I mentioned as an exemple fluid simulation. The intern then asked for the repository and I had nothing to show for it, so I decided it was time to start over.
I wanted to have the simulation in Fortran, but since this is the kind of code that is extremely hard to debug, I planned on doing it in Octave first.
I reinstalled Octave on my compurer and within a few hours I was aggravated. I cannot believe I was able to work with such limitations in the past. Octave (and MATLAB) can have only one function per file. There is basically no possibility to have any types specified and everything is backwards. I decided to port my code, before it was too late, to something less appalling. I thought about Python, but I had a flash that Julia could support something like that and that it also had a built-in linear solver and good visualisation tools. I also wanted to learn the language. So I installed Julia, ported what little code I had already, and was almost instantly rewarded. The only thing that is annoying is that the visualisation tools are not as good as MATLAB's, the libraries are not as mature. But at least it is workable. More information on Julia.
Bridson's book covers this part in Chapter 5. Chapter 4 is for collisions and boundaries, so I read it, but did not implement anything from it.
I won't go into mathematical details, but the book ends up with an equation of the form \( \boldsymbol{A} \cdot \vec{p} = \vec{b} \), the classic form of a linear system. Now I feel like a physicist again! The vector \(\vec{p}\) is the pressure at each cell middle in the grid, so it is of length \(N^2\) if the grid is a square with a width of N cells. This means that the \(\vec{b}\) vector is of the same length, and that the matrix \( \boldsymbol{A} \) is a huge matrix of size \( N^4 \). So a 10 by 10 grid will generate 100 pressure values, but a matrix with 10'000 values. A 100 by 100 grid, which is reasonable, would yield a matrix with 100'000'000 values!
I made my first tests with full matrices, but used built-in sparse matrices from Julia soon after. In the future, if I want custom algorithms and sparse matrices in Fortran, I'll probably need to support them myself. I was able to build the matrix using section 5.3 of the book and put equations found on page 74 into a matrix and a vector. The vector \(\vec{b}\) depends on velocities in the grid, but the matrix \( \boldsymbol{A} \) itself only depends on a few constants and the number of neighbors a cell has, so it can be constructed only once for the whole duration of the simulation and reused each frame.
The good thing with Julia, is that I can build the matrix and vector, and then just ask it to solve the problem without coding anything, just like in MATLAB. It is quite slow, but will yield a valid result. In the following example, I put velocity only in X and going only on the right as inputs for the solver, which should not be allowed by physics.
The solver should fix the issue and find a physical solution where the pressure is related to velocity, but also where there is no air (I simulate air) that is created, so no air can enter or exit the simulation. The pressure result is given on next picture.
The last step of the solver is actually to apply the pressure as a force to modify velocity. This is how we will know if the result actually makes sense. Next picture gives it, and it looks like something. It now has Y velocity, and modified X velocity, but aside from looking very smooth it's hard to understand what is happening still.
I decided to show the same result, but using arrows to represent velocities and have an actual understanding of what is happening. I think that in the next picture, it is clearer that the input velocity it set does not make sense. It appears out of nowhere, and stops out of nowhere in a disk. Note that there is no guarantee that the solver can actually solve a system that is very wrong from the start, and a valid solution should be provided as a starting point.
The result of the simulation, given in terms of arrows actually shows what was hard to visualize with the previous graphs. The solver actually tried to keep the X velocity as much as possible in the center, but has to loop the velocity arround so no air is created, but actually flows in circular movements. If air was entering the system by the boundary, the result would be different, but here nothing enters and nothing exits.
Now that I proved that the simulation was possible, and quite fast in Julia, the solver needs to be specialized for performance. The book helps us with doing just that, giving us a very specialized and extremely fast algorithm that I did not have time to implement yet. Instead, I tried to do a somewhat generic solver using the Gaussian Elimination algorithm. This is the algorithm that we all learned in college, for more information : Gaussian elimination.
Eventually I will have to implement the real algorithm, but I wanted to experiment a bit with the solver. It is useful, because it allows me to also implement sparse matrix custom implementation. There are multiple ways to handle sparse matrices. The simplest is just to have a list of structures, each containing a column, a row and a value. In 2D, my algorithm uses a maximum of 5 values per row of the matrix, and each row contains one value for each cell of the grid. So if our grid contains 10'000 cells, our matrix is at a maximum filled at a rate of 5/10'000, so 0.05%. Even by storing both the row and the column we save a huge amount of memory. I implemented the basic operations on this type of matrix and wrote some unit tests for it (like multiplications). I also made it work with the solver, but the way I insert values makes it very slow.
Another way to store sparse matrix is to use the compressed sparse column structure (or compressed sparse row) where values of the same column (or row) are all put together in memory, so we only need to keep track of how many there are in that row or column, saving even more in memory, and making the memory access faster. It is explained here . I actually had to fill manually this type of matrix (compressed sparse column, because MATLAB works that way) in the past, so implementing my own (compressed sparse row I think) will be easy enough.
I am now in the process of converting the solver to these types of matrices, but I am still unfinished. After I am done and have working custom sparse matrices, I will implement Bridson's algorithm. For now, it takes about 33ms in Julia on one thread using my bad solver for a 40x40 grid. It looks like it is quite viable.
Starting from now, I'll add the stats of the codebase for the engine.
-------------------------------------------------- ['cpp'] 56 files and 6852 lines ['h', 'inl'] 75 files and 4271 lines ['cs'] 11 files and 2049 lines ['rs'] 10 files and 1562 lines ['jl'] 14 files and 1350 lines ['zig'] 10 files and 562 lines ['glsl'] 9 files and 189 lines ['py'] 3 files and 188 lines ['ps1', 'sh'] 10 files and 142 lines ['cmake', 'CMakeLists.txt'] 5 files and 126 lines ['ispc'] 1 files and 18 lines ['f90'] 1 files and 13 lines -------------------------------------------------- Total : 205 files and 17322 lines --------------------------------------------------